//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/ephi/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>

#include "ephi.hpp"

static inline bool pis (prec_t a, prec_t b) { return a == b || fabs(a - b) < 1e-10; }

int main (int argc, char *argv[])
{
  {
  LinearSystem ls(3);
  // x+y+z=3
  ls.getRow(0)[0] = 1;
  ls.getRow(0)[1] = 1;
  ls.getRow(0)[2] = 1;
  ls.getRow(0)[3] = 3;
  // 2x+3y+7z=0
  ls.getRow(1)[0] = 2;
  ls.getRow(1)[1] = 3;
  ls.getRow(1)[2] = 7;
  ls.getRow(1)[3] = 0;
  // x+3y-2z=17
  ls.getRow(2)[0] = 1;
  ls.getRow(2)[1] = 3;
  ls.getRow(2)[2] = -2;
  ls.getRow(2)[3] = 17;

  if ( !ls.solveGBS() )
    {
      printf ("NOK: solving solution one failed!\n");
    }
  else
    {
      prec_t x = ls.getResults()[0];
      prec_t y = ls.getResults()[1];
      prec_t z = ls.getResults()[2];

      if ( pis(x, 1) && pis(y, 4) && pis(z, -2) )
          printf ("OK: GBS results are: x=%.16f, y=%.16f, z=%.16f\n", prec2double(x), prec2double(y), prec2double(z));
      else
          printf ("NOK: GBS wrong results: x=%.16f, y=%.16f, z=%.16f\n", prec2double(x), prec2double(y), prec2double(z));
    }
  }
  {
  LinearSystem ls(3);
  // x+y+z=3
  ls.getRow(0)[0] = 1;
  ls.getRow(0)[1] = 1;
  ls.getRow(0)[2] = 1;
  ls.getRow(0)[3] = 3;
  // 2x+3y+7z=0
  ls.getRow(1)[0] = 2;
  ls.getRow(1)[1] = 3;
  ls.getRow(1)[2] = 7;
  ls.getRow(1)[3] = 0;
  // x+3y-2z=17
  ls.getRow(2)[0] = 1;
  ls.getRow(2)[1] = 3;
  ls.getRow(2)[2] = -2;
  ls.getRow(2)[3] = 17;

  if ( !ls.solveGJ() )
    {
      printf ("NOK: solving solution one failed!\n");
    }
  else
    {
      prec_t x = ls.getResults()[0];
      prec_t y = ls.getResults()[1];
      prec_t z = ls.getResults()[2];

      if ( pis(x, 1) && pis(y, 4) && pis(z, -2) )
          printf ("OK: GJ results are: x=%.16f, y=%.16f, z=%.16f\n", prec2double(x), prec2double(y), prec2double(z));
      else
          printf ("NOK: GJ wrong results: x=%.16f, y=%.16f, z=%.16f\n", prec2double(x), prec2double(y), prec2double(z));
    }
  }
  {
  LinearSystem ls(3);
  // x+y+z=1
  ls.getRow(0)[0] = 1;
  ls.getRow(0)[1] = 1;
  ls.getRow(0)[2] = 1;
  ls.getRow(0)[3] = 1;
  // 2x+y+z=3
  ls.getRow(1)[0] = 2;
  ls.getRow(1)[1] = 1;
  ls.getRow(1)[2] = 1;
  ls.getRow(1)[3] = 3;
  // x+y-2z=4
  ls.getRow(2)[0] = 1;
  ls.getRow(2)[1] = 1;
  ls.getRow(2)[2] = -2;
  ls.getRow(2)[3] = 4;

  if ( !ls.solveGBS() )
    {
      printf ("NOK: solving solution one failed!\n");
    }
  else
    {
      prec_t x = ls.getResults()[0];
      prec_t y = ls.getResults()[1];
      prec_t z = ls.getResults()[2];

      if ( pis(x, 2) && pis(y, 0) && pis(z, -1) )
          printf ("OK: GBS results are: x=%.16f, y=%.16f, z=%.16f\n", prec2double(x), prec2double(y), prec2double(z));
      else
          printf ("NOK: GBS wrong results: x=%.16f, y=%.16f, z=%.16f\n", prec2double(x), prec2double(y), prec2double(z));
    }
  }
  {
  LinearSystem ls(3);
  // x+y+z=1
  ls.getRow(0)[0] = 1;
  ls.getRow(0)[1] = 1;
  ls.getRow(0)[2] = 1;
  ls.getRow(0)[3] = 1;
  // 2x+y+z=3
  ls.getRow(1)[0] = 2;
  ls.getRow(1)[1] = 1;
  ls.getRow(1)[2] = 1;
  ls.getRow(1)[3] = 3;
  // x+y-2z=4
  ls.getRow(2)[0] = 1;
  ls.getRow(2)[1] = 1;
  ls.getRow(2)[2] = -2;
  ls.getRow(2)[3] = 4;

  if ( !ls.solveGJ() )
    {
      printf ("NOK: solving solution one failed!\n");
    }
  else
    {
      prec_t x = ls.getResults()[0];
      prec_t y = ls.getResults()[1];
      prec_t z = ls.getResults()[2];

      if ( pis(x, 2) && pis(y, 0) && pis(z, -1) )
          printf ("OK: GJ results are: x=%.16f, y=%.16f, z=%.16f\n", prec2double(x), prec2double(y), prec2double(z));
      else
          printf ("NOK: GJ wrong results: x=%.16f, y=%.16f, z=%.16f\n", prec2double(x), prec2double(y), prec2double(z));
    }
  }
  return 0;
}

